In [97]:
import pandas
from datetime import datetime, timedelta
import matplotlib.pyplot as plt
In [9]:
%matplotlib inline
pandas.set_option("display.max_rows", 10)
SunPy Figure
In [10]:
import sunpy.map
import sunpy
aiamap = sunpy.map.Map(sunpy.AIA_171_IMAGE)
smap = aiamap.submap([-1200, -200], [-1000, -0])
from sunpy import lightcurve
from sunpy.time import TimeRange
goes = lightcurve.GOESLightCurve.create(TimeRange('2012-07-12 14:00', '2012-07-12 23:00'))
compmap = sunpy.map.Map(sunpy.AIA_171_IMAGE, sunpy.RHESSI_IMAGE, composite=True)
compmap.set_levels(1, range(0, 50, 5), percent=True)
compmap.set_colors(1, 'Reds_r')
fig = plt.subplots(ncols=2)
plt.subplot(231)
smap.plot()
smap.draw_grid()
smap.draw_limb()
ax = plt.subplot(232)
compmap.plot()
ax.axis([200, 600, -600, -200])
plt.subplot(233)
ax1 = goes.plot()
ax1.set_yscale('log')
plt.subplots_adjust(wspace=0.5)
plt.show()
In [4]:
from sunpy.net import hek
client = hek.HEKClient()
In [5]:
tstart = '2012/07/12 07:23:56'
tend = '2012/07/14 12:40:29'
event_type = hek.attrs.EventType('FL')
conditions = (hek.attrs.OBS.Instrument == 'GOES') & (hek.attrs.FRM.Name == 'SWPC') & (hek.attrs.FL.GOESCls.like('X%') | hek.attrs.FL.GOESCls.like('M%'))
result = client.query(hek.attrs.Time(tstart,tend), event_type, conditions)
In [6]:
len(result)
Out[6]:
In [7]:
for res in result:
print(res.get('fl_goescls'))
In [8]:
hek_flare = result[0]
In [9]:
hek_flare.get('fl_goescls')
Out[9]:
In [10]:
hek_flare.get('hpc_x'), hek_flare.get('hpc_y')
Out[10]:
In [11]:
tstart = '2012/07/12 07:23:56'
tend = '2012/07/15 12:40:29'
event_type = hek.attrs.EventType('CE')
conditions = conditions = hek.attrs.OBS.Instrument == 'SOHO'
result = client.query(hek.attrs.Time(tstart,tend), event_type)
In [12]:
result
Out[12]:
In [13]:
from sunpy.net import hek2vso
Set up some nice plotting for time series data
In [14]:
import seaborn as sns
sns.set(style="darkgrid")
In [15]:
file = '/Users/schriste/Desktop/SunPy_ODDE_Proposal/AC_H1_EPM_201207.txt'
In [16]:
ace_col_names = ['FP6P_.761-1.22MEV_IONS', 'FP7P_1.22-4.97MEV_IONS', 'UNC_FP6P_.761-1.22MEV_IONS', 'UNC_FP7P_1.22-4.97MEV_IONS']
names = ['date', 'hour'] + ace_col_names
ace_data = pandas.read_csv(file, skiprows=74, delim_whitespace=True, names = names)
In [17]:
ace_data = ace_data.truncate(after=len(ace_data)-5)
In [18]:
times = [datetime.strptime(t[0:-4], '%d-%m-%Y %H:%M:%S') for t in ace_data['date'] + ' ' + ace_data['hour']]
In [19]:
ace_data['times'] = times
ace_data = ace_data.set_index('times')
In [20]:
ace_data = ace_data.drop('hour',1)
ace_data = ace_data.drop('date',1)
In [21]:
for col in ace_col_names:
ace_data[col] = ace_data[col].convert_objects(convert_numeric=True)
In [22]:
ace_data.dtypes
Out[22]:
In [23]:
for col in ace_col_names:
ace_data[col][ace_data[col] < 0] = np.nan
In [24]:
palette = sns.color_palette()
In [25]:
plt.figure()
ace_data[ace_col_names[0]].plot()
ace_data[ace_col_names[1]].plot()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.title("EPAM/ACE Electron Proton Alpha Monitor")
plt.ylabel("Ion Flux")
plt.show()
In [26]:
def read_iswa_enlil_dat(file, col_names):
names = ['year', 'month', 'day', 'hour', 'minute'] + col_names
data = pandas.read_csv(file, skiprows=2, names=names, delim_whitespace=True)
times_str = [str(int(t[1]['year'])) + '-' + str(int(t[1]['month'])).zfill(2) + '-' + str(int(t[1]['day'])).zfill(2) + ' ' + str(int(t[1]['hour'])).zfill(2) + ':' + str(int(t[1]['minute'])).zfill(2) for t in data.iterrows()]
times = [datetime.strptime(t, '%Y-%m-%d %H:%M') for t in times_str]
data['times'] = times
data = data.set_index('times')
data = data.drop('year',1)
data = data.drop('month',1)
data = data.drop('day',1)
data = data.drop('hour',1)
data = data.drop('minute',1)
return data
In [27]:
def read_iswa_data_dat(file):
data = pandas.read_csv(file, delim_whitespace=True)
times_str = [str(t[0]) + ' ' + str(t[1]['Timestamp'])[0:-2] for t in data.iterrows()]
times = [datetime.strptime(str_time, '%Y-%m-%d %H:%M:%S') for str_time in times_str]
data['times'] = times
data = data.set_index('times')
data = data.drop('Timestamp',1)
return data
In [28]:
file = '/Users/schriste/Desktop/20120712_193500_ENLIL_time_line.dat.txt'
In [29]:
col_names = ['B_enl','V_enl','n_enl','T_enl']
enlil_data = read_iswa_enlil_dat(file, col_names=col_names)
In [30]:
enlil_data
Out[30]:
In [31]:
plt.figure()
enlil_data.plot(subplots=True, figsize=(10, 10))
plt.show()
Deriving the Kp index from ENLIL data
In [32]:
kp90exp=9.5-np.exp(2.17676-(0.000052001)*((enlil_data.V_enl)**1.3333)*((enlil_data.B_enl)**0.6666)*(np.sin((np.pi/2)/2)**(8/3.)))
kp135exp=9.5-np.exp(2.17676-(0.000052001)*((enlil_data.V_enl)**1.3333)*((enlil_data.B_enl)**0.6666)*(np.sin((3*np.pi/4)/2)**(8/3.)))
kp180exp=9.5-np.exp(2.17676-(0.000052001)*((enlil_data.V_enl)**1.3333)*((enlil_data.B_enl)**0.6666)*(np.sin((np.pi)/2)**(8/3.)))
In [33]:
plt.figure()
kp90exp.plot()
kp135exp.plot()
kp180exp.plot()
plt.title("ENLIL Predicted KP Index")
plt.ylabel("KP")
plt.show()
In [34]:
file = '/Users/schriste/Desktop/SunPy_ODDE_Proposal/20120712_193500_ENLIL_Kp_timeline.dat.txt'
In [35]:
kp_col_names = ['Kp_90', 'Kp_135', 'Kp_180']
enlil_kp_data = read_iswa_enlil_dat(file, col_names=kp_col_names)
In [36]:
enlil_kp_data
Out[36]:
In [37]:
plt.figure()
for col in kp_col_names:
enlil_kp_data[col].plot()
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.title("ENLIL Predicted KP Index")
plt.ylabel("KP")
plt.show()
a current bug in our GOES data fetcher is that it cannot handle multiple days in the time range
In [38]:
from sunpy.lightcurve import GOESLightCurve
from sunpy.time import TimeRange
In [39]:
tr = TimeRange(ace_data.index[0].to_datetime(), ace_data.index[-1].to_datetime())
In [40]:
tr
Out[40]:
In [41]:
goes = GOESLightCurve.create(tr)
In [42]:
plt.figure()
ax = goes.data['xrsa'].plot()
plt.yscale('log')
ax = goes.data['xrsb'].plot()
plt.show()
note to self: had to manually delete blank lines at the end of file otherwise times have nan's and fails to parse
In [45]:
file = '/Users/schriste/Desktop/SunPy_ODDE_Proposal/ISWA_GOES.txt'
In [46]:
iswa_goes = read_iswa_data_dat(file)
In [47]:
iswa_goes
Out[47]:
In [48]:
iswa_goes['Long_Wave'][iswa_goes['Long_Wave'] < 0] = np.nan
In [49]:
plt.figure()
iswa_goes.plot()
plt.yscale('log')
plt.title('GOES 1-8 A (ISWA)')
plt.ylabel('Watts')
plt.show()
In [52]:
file = '/Users/schriste/Desktop/SunPy_ODDE_Proposal/ISWA_KP.txt'
In [53]:
iswa_kp = read_iswa_data_dat(file)
In [54]:
plt.figure()
iswa_kp.plot()
plt.title('KP (ISWA)')
plt.ylabel('KP Index')
plt.show()
In [55]:
type(iswa_kp['KP'])
Out[55]:
In [56]:
plt.figure()
iswa_kp.plot()
kp135exp.plot()
kp180exp.plot()
kp90exp.plot()
plt.legend()
plt.show()
In [137]:
kpiswa = iswa_kp['KP'].reindex(enlil_kp_data['Kp_90'].index, method='ffill')
kp = pandas.DataFrame({"ISWA":kpiswa, "ENLIL":enlil_kp_data['Kp_90'], "Obs Mean": pandas.rolling_mean(kp['ISWA'], 120)})
In [138]:
kp
Out[138]:
In [139]:
kp.corr()
Out[139]:
In [141]:
plt.figure()
ax = kp.plot()
ax.ylabel("K-index")
plt.show()
In [136]:
kp.corr()
Out[136]:
In [61]:
plt.figure()
plt.plot(kp['ISWA'], kp['ENLIL'], '*')
plt.show()
Documentation for the following file http://www.srl.caltech.edu/STEREO/Public/HET_public.html
In [65]:
file = '/Users/schriste/Desktop/SunPy_ODDE_Proposal/BeH12Jul.1m.txt'
In [66]:
def read_stereob_impact(file):
names = ['0', 'year', 'month', 'day', 'hhmm',
'eflux_0.7-1.4', 'eflux_0.7-1.4_uncert',
'eflux_1.4-2.8', 'eflux_1.4-2.8_uncert',
'eflux_2.8-4.0', 'eflux_2.8-4.0_uncert',
'pflux_13.6-15.1', 'pflux_13.6-15.1_uncert',
'pflux_14.9-17.1', 'pflux_14.9-17.1_uncert',
'pflux_17.0-19.3', 'pflux_17.0-19.3_uncert',
'pflux_20.8-23.8', 'pflux_20.8-23.8_uncert',
'pflux_23.8-26.4', 'pflux_23.8-26.4_uncert',
'pflux_26.3-29.7', 'pflux_26.3-29.7_uncert',
'pflux_29.5-33.4', 'pflux_29.5-33.4_uncert',
'pflux_33.4-35.8', 'pflux_33.4-35.8_uncert',
'pflux_35.5-40.5', 'pflux_35.5-40.5_uncert',
'pflux_40.0-60.0', 'pflux_40.0-60.0_uncert',
'pflux_60.0-100.0', 'pflux_60.0-100.0_uncert']
data = pandas.read_csv(file, delim_whitespace=True, skiprows=24, dtype={'hhmm': np.dtype('str')}, header=None, names=names[0:40])
dates_str = [str(row[1][1]) + '-' + str(row[1][2]).zfill(1) + '-' + str(row[1][3]).zfill(2) + ' ' + row[1][4][0:2] + ':' + row[1][4][2:4] for row in data.iterrows()]
times = [datetime.strptime(str_time, '%Y-%b-%d %H:%M') for str_time in dates_str]
data = data.drop('year', 1)
data = data.drop('month', 1)
data = data.drop('hhmm', 1)
data = data.drop('day', 1)
data = data.drop('0', 1)
data['times'] = times
data = data.set_index('times')
return data
In [67]:
data = read_stereob_impact(file)
In [68]:
cols = data.columns
In [69]:
data
Out[69]:
In [70]:
e_cols = ['eflux_0.7-1.4', 'eflux_1.4-2.8', 'eflux_2.8-4.0']
p_cols = ['pflux_13.6-15.1', 'pflux_14.9-17.1', 'pflux_17.0-19.3', 'pflux_20.8-23.8', 'pflux_23.8-26.4',
'pflux_26.3-29.7', 'pflux_29.5-33.4', 'pflux_33.4-35.8', 'pflux_35.5-40.5', 'pflux_40.0-60.0',
'pflux_60.0-100.0', ]
e_uncert_cols = [e_col + '_uncert' for e_col in e_cols]
p_uncert_cols = [p_col + '_uncert' for p_col in p_cols]
In [71]:
e_label = ['Electron Flux ' + e_col.split('_')[1] + ' keV' for e_col in e_cols]
p_label = ['Proton Flux ' + p_col.split('_')[1] + ' keV' for p_col in p_cols]
In [72]:
plt.figure()
for col, label in zip(e_cols, e_label):
data[col].plot(label=label)
plt.legend(loc=2)
plt.title('Electrons')
plt.show()
plt.figure()
for col, label in zip(p_cols, p_label):
data[col].plot(label=label)
plt.title('Protons')
plt.legend(loc=2)
plt.show()
In [73]:
data.resample('T')
plt.figure()
for col in e_cols:
data[col].resample('H', how='mean').plot()
plt.legend()
plt.title('Electrons')
plt.show()
In [74]:
time_range = ['2012-07-12 00:00:00', '2012-07-15 00:00:00']
col_index = 2
col = p_cols[col_index]
col_uncert = p_uncert_cols[col_index]
label = p_label[col_index]
p_zoom = data[col].truncate(before=time_range[0], after=time_range[1])
p_uncert_zoom = data[col_uncert].truncate(before=time_range[0], after=time_range[1])
p_plus = p_zoom + p_uncert_zoom
p_minus = p_zoom - p_uncert_zoom
p_plus = p_plus.resample('H', how='max')
p_minus = p_minus.resample('H', how='min')
In [91]:
plt.rc('text', usetex=True)
plt.figure()
p_zoom.resample('H', how='mean').plot(style='k', label=label)
plt.fill_between(p_plus.index, p_plus.values, y2=p_minus.values, interpolate=True, alpha=0.5)
plt.title("STEREO B IMPACT")
plt.ylabel(r'particles cm$^{-2}$ sr$^{-1}$ sec$^{-1}$ MeV$^{-1}$')
plt.legend()
plt.show()
warning! the widget below will only work in live notebooks
In [3]:
import sunpy.map
from IPython.html.widgets import interact, interactive
import matplotlib.pyplot as plt
from IPython.display import display
%matplotlib inline
In [4]:
aia = sunpy.map.Map(sunpy.AIA_171_IMAGE)
In [5]:
def plot_map(g=3, vmin=0, vmax=1000):
fig = plt.figure()
aia.cmap.set_gamma(g)
aia.mpl_color_normalizer.vmin=vmin
aia.mpl_color_normalizer.vmax=vmax
aia.plot()
plt.colorbar()
plt.show()
The following is an interactive widget which lets you set the gamma of the plot interactively.
In [6]:
w = interactive(plot_map, g=(0.0,1.0,0.1), vmin=(0.0, aia.max(), 10), vmax=(0.0, aia.max(), 10))
display(w)
This is a bit slow. Let's try to speed it up but plotting less data (i.e. smaller map)
In [7]:
aia = sunpy.map.Map(sunpy.AIA_171_IMAGE).submap([-1200,-900], [-400, -200])
aia
Out[7]:
In [8]:
w = interactive(plot_map, g=(0.0,5.0,0.1), vmin=(0.0, aia.max(), 10), vmax=(0.0, aia.max(), 10))
display(w)
By messing with the sliders, we can easily pick out a faint loop in the corona
In [93]:
file = '/Users/schriste/Desktop/SunPy_ODDE_Proposal/pointdata_27423.txt'
In [112]:
def read_ccmc(file):
names = ['days', 'R', 'Lat', 'Lon', 'V_r', 'V_lon', 'V_lat', 'B_r', 'B_lon', 'B_lat', 'N', 'T', 'E_r', 'E_lon', 'E_lat', 'V', 'B', 'P_ram', 'BP']
data = pandas.read_csv(file, delim_whitespace=True, skiprows=7, header=None, names=names)
tstart = datetime.strptime('2002/03/02 00:00:00', '%Y/%m/%d %H:%M:%S')
time = [tstart + timedelta(row[1]['days'],0,0) for row in data.iterrows()]
data.index = time
data = data.drop('days', 1)
return data
In [115]:
ccmc_enlil = read_ccmc(file)
In [134]:
plt.figure(1)
plt.subplot(311)
ax = ccmc_enlil['V'].plot()
ax.set_xticklabels([])
ax.set_title('CCMC ENLIL')
plt.ylabel(r'V [km s$^{-1}$')
plt.subplot(312)
ax = ccmc_enlil['B'].plot()
plt.ylabel(r'B [nT]')
ax.set_xticklabels([])
plt.subplot(313)
ax = ccmc_enlil['N'].plot()
plt.ylabel(r'N [cm$^{-3}$]')
plt.show()
In [ ]: